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■ A set of energy-momentum constraint-preserving boundary conditions is proposed for the 
' first-order Z4 case. The stability of a simple numerical implementation is tested in the linear 

fN| \ regime (robust stability test), both with the standard corner and vertex treatment and with a 

^P, modified finite-differences stencil for boundary points which avoids corners and vertices even 

I in cartesian-like grids. Moreover, the proposed boundary conditions are tested in a strong 

\^ • field scenario, the Gowdy waves metric, showing the expected rate of convergence. The 

I— i' accumulated amount of energy- momentum constraint violations is similar or even smaller 

than the one generated by either periodic or refiection conditions, which are exact in the 
Gowdy waves case. As a side theoretical result, a new symmetrizer is explicitly given, 
which extends the parametric domain of symmetric hyperbolicity for the Z4 formalism. The 
I application of these results to first-order BSSN-like formalisms is also considered. 
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2 ! I. INTRODUCTION 

o : 

■ Constraint-preserving boundary conditions is an active research topic in Numerical Relativity. 
^ , During the first half of this decade, many conditions have been proposed, adapted in each case to 

c5 ■ some specific 3-M evolution formalism: Fritelli-Reula |J], KST [3, B], BSSN-NOR [3], or Z4 Q]. 

The focus changed suddenly after 2005 by the impact of a breakthrough: the first ' long term' 
binary-black-hole simulation, which was achieved in a generalized-harmonic formalism [sjl. A series 
of constraint-preserving boundary conditions proposals in this framework started then 7, ^, and 
continues today [l^ . 

We will retake in this paper the 3-1-1 approach to constraint-preserving boundary conditions, 



following the way opened very recently for the BSSN case ll|]. More specifically, we will revisit 



the Z4 case, not just because of its intrinsic relevance, but also for its relationship with other 3-1-1 
formulations (BSSN, KST, see refs. {l3. \\\ for details). Also, the close relationship between the 
Z4 and the generalized-harmonic formulations suggest that our results could provide a different 
perspective in this other context. This was actually what happened with the current constraint- 



damping terms: first derived in the Z4 context [14'] and then appHed successfully in generalized- 
harmonic simulations 

Our results are both at the theoretical and the numerical level. In section II, we consider the 
first-order Z4 formalism in normal coordinates (zero shift) for the harmonic slicing case. This case 
was known to be symmetric- hyperbolic for a particular choice of the parameter which controls the 
ordering of space derivatives [la. [l6|. We extend this result to a range of this ordering parameter, 
by providing explicitly a positive-definite energy estimate. Then we use this estimate for deriving 
algebraic constraint-preserving boundary conditions both for the energy and the normal momentum 
components. 

In section III we consider the dynamical evolution of constraint violations (subsidiary system). 
Following standard methods [2-4], we transform algebraic boundary conditions of the subsidiary 
system into derivative boundary conditions for the main system. We introduce a new basis of 
dynamical fields in order to revise the constraint-preserving conditions proposed in refs. P, [3] for 
the Z4 formalism, including also a new coupling parameter which affects the propagation speeds of 
the (modified) incoming modes. In the case of the energy constraint, we get a closed subsystem for 
the principal part, allowing an analytical stability study at the continuum level which is presented 
in Appendix B. 

A simple numerical implementation of the proposed conditions is given in section IV, where we 
test the stability in the linear regime, by considering small random- noise perturbations aroun flat 
space (robust stability test). The results show the numerical stability of the proposed boundary 
conditions in this regime for many different combinations of the parameters. The space discretiza- 
tion scheme is the simplest one with the summation-by-parts (SBP) property [17]. In this way 
we avoid masking the effect of our conditions (at the continuum level) with the effect of more ad- 



vanced space-discretization algorithms, like FDOC ^8|] devised to reduce the high frequency noise 
level in long-term simulations, which has recently been applied to the black-hole case [3]. For a 
comparison, we run also with periodic boundary conditions, where the noise level keeps constant. 
The proposed boundary conditions produce instead a very effective decreasing of (the cumulated 
effect of) energy and momentum constraint violations. In the case of cartesian-like grids, we also 
compare the standard 'a la Olsson' treatment [17], with a modified numerical implementation 
which does not use the corner and vertex points, avoiding in this way some stability issues and 
providing much cleaner evidence of constraint preservation. 

In section V, we test the non-linear regime with the Gowdy waves [2^ metric, one of the standard 
numerical relativity code tests, as we have done recently for the energy constraint case 2]|. We 
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endorse in this w ay s ome recent claims (by Winicour and others) that the current code cross- 



comparison efforts 



22 



23] should be extended to the boundaries treatment. A convergence test is 



performed against this exact strong-field solution, showing the expected convergence rate (second 
order for our simple SBP method). Testing the proposed boundary conditions results into a stable 
and constraint-preserving behavior, in the sense that energy and momentum constraint violations 
remain similar or even smaller than the corresponding effects with exact (periodic or reflection) 
boundary conditions for the Gowdy metric. 



II. THE Z4 CASE REVISITED 



We will consider here the Z4 evolution system: 



(1) 



More specifically, we will consider the first-order version in normal coordinates, as described 



in refs. 
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la ]. For further convenience, we will recombine the basic first-order fields 



{Kij, Dijk, Ai, 6, Zi) in the following way: 



Vi 
Wi 



Kij -{trK- 9) 

Y^i^irs — Dris) — Zk 
Dijk - (Y^Dirs - Vi)-fjk 
Ai-Y'Dirs + 2Vi 



(2) 
(3) 
(4) 
(5) 



so that the new basis is {Hij, iJ-ijk, Wi, Q, Vi). Note that the vector Zi can be recovered easily as 



ik- 



(6) 



With this new choice of basic dynamical fields, the principal part of the evolution system gets 
a very simple form in the harmonic slicing case: 



dtWi 
dtQ + adkV'' 
dtVi + adiQ 
dt Uij + a dk X^ij 
dt fikij + a dk Uij 



(7) 
(8) 
(9) 
(10) 

(11) 
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where the dots stand for non-principal contributions, and we have noted for short 

Afcjj = ^ikij + ik(iWj) - Wk jij - (1 + C) [ fj-{ij)k + ikiiZj) ] , (12) 

where C is a space-derivatives ordering parameter and round brackets denote index symmetriza- 
tion. 



The first-order version of the Z4 system is known to be symmetric- 
dinates with harmonic shcing, at least for the usual ordering ^ = — 1 



lyperbolic in normal coor- 



]J]. It follows from ifTlfTT]) 



that this result can be extended to the following range of the ordering parameter 

-1<C<0, (13) 

which covers the symmetric ordering case = 0). The corresponding 'symmetrizer', or 'energy 
estimate', can be written as: 

s = e^ + VkV'^ + u'^ + fi^'^jikij + (1 + c) {z'^Zk - jJ^'^hjk) + 2CZkW\ (i4) 

where we have noted for short 

Jj'kij = fJ-kij - Wk 'Jij ■ (15) 

Allowing for ([TlfTTI). we get 

^dtS = dkiQv'' + u,jX^'^) + --- (16) 

and the divergence theorem can be used in order to complete the proof. The positivity proof for 
S for the interval (I13p is given in Appendix A. 



Characteristic decomposition 

We can consider now some specific space surface, in order to identify the constraint modes by 
looking at the evolution equations for B and Zi in the system (fTl llip . It follows from ([HI [9]) that 
the energy-constraint modes are given by the pair 

E^ = e ± Vn (17) 

with propagation speed ±a (the index n meaning the projection along the unit normal rii). Also, 
allowing for (I6|lip . we can easily recover the evolution equation for Zi, namely 

dtZi-adkU''i = --- (18) 
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so that we can identify the momentum-constraint modes with the three pairs, with propagation 
speed ±a, 

M± = Uni ± Xnni ■ (19) 

Note that, allowing for (l2|), the normal component linn does correspond with the transverse-trace 
component of the extrinsic curvature Kij. We give for completeness the remaining modes, the fully 
tangent ones, with propagation speed ±a, 

^AB = ^AB ± XnAB (20) 

(capital indices denote a projection tangent to the surface), and the standing modes (zero propa- 
gation speed): 

, Va, ^iA^j ■ (21) 



Algebraic boundary conditions 

We can take advantage of the positive-definite energy estimate (jl4p in order to derive suitable 
algebraic boundary conditions. We can integrate (jl6p in space and, by applying the divergence 
theorem, we get a positivity condition for the boundary terms, namely 

(rf^' Xnij + e K) Is > (22) 

where S stands for the boundary surface (n being here its outward normal). 

The contribution of the fully tangent modes (I20p . independent of the energy and momentum 
sectors, is given by 

n^"" XnAB = \tr [{T+ f - [T-f] , (23) 

so that the contribution of these modes to the boundary term in (I22p will be non-negative if we 
impose the standard algebraic boundary-conditions: 

Tab = '^T+^ kl<l' (24) 

the case a = corresponding to maximal dissipation. A less strict condition is obtained by adding 
an inhomogeneous term, namely 

TXb = <'T+^ + Gab. (25) 
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This can cause some growth of the energy estimate but, provided that the array G consists of 
prescribed spacetime functions, the growth rate can be bounded in a suitable way so that a well- 



posed system can still be obtained (see for instance refs. |4l. 

This simple strategy, when applied to the energy and momentum modes (|17l [T9|) is not com- 
patible with constraint preservation in the generic case (see also ref. [3]). For the energy sector, 
constraint preservation is obtained only for the extreme case: 

= E+ ^ G = , (26) 

which will reflect energy-constraint violations back into the evolution domain. These conditions 
would be then of a limited practical use in realistic simulations. 

A different approach can be obtained by realizing that the contribution to the boundary term 
in (j22]) would have the right sign if one uses the following ' logical gate' condition: 

e Is = if (9 K) Is < (27) 

(0-gate in ref. jsij]). It is clear that the boundary condition (j27p preserves the energy constraint, 
as it modifies just the Q values, by setting them to zero when the condition is fulfilled, without 
affecting any other dynamical field. 

The same strategy can work for normal components of the momentum modes (jl9p . at least for 
the symmetric choice of the ordering parameter. Allowing for (|12p . one has 

M± = n„„ ^Zn (C = 0) , (28) 

so that a constraint-preserving (reflection) condition can be obtained in the extreme case as well. 
In the logical gate approach, the contribution of the modes ()28p to the boundary term in ()22p will 
have the right sign if one uses the condition (case C = only): 

Zn Is =0 if {Zn Unn) |e > , (29) 

which clearly preserves the normal component of the momentum constraint. 

For the tangent momentum modes (tangent to the boundary surface), however, the contri- 
bution in ([22]) win be 

2n"^A„n^, (30) 

where XnnA is inhomogeneous in Za for any value of the ordering parameter. Moreover, the 
inhomogeneous terms are not prescribed functions, but rather some combinations of dynamical 
fields. A different strategy must then be devised in this case, as we will see below. 
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III. CONSTRAINTS EVOLUTION AND DERIVATIVE BOUNDARY CONDITIONS 

The time evolution of the energy-momentum constraints can be easily derived by taking the 
divergence of the Z4 field equations ([T]) , that is 

UZ^ + R^.Z'' = . (31) 

We can write down the second order equation ()3ip as a first order system and impose then maxi- 
mally dissipative boundary conditions on (the first derivatives of) the components. In this way, 
the boundaries will behave as one-way membranes for constraint-violating modes, at least for the 
ones propagating along the normal direction nj. 

The procedure can be illustrated with the energy-constraint, that is the time component of ()3ip : 



4 G - A G = • • • (32) 

A first-order version can be obtained as usual by considering first-order derivatives as independent 
quantities, namely 

G = 1/a a* Q , Qk = dkQ. (33) 

We can write then (|32p as the following first-order symmetric-hyperbolic system 

dtQ-adkQ'' = ■■■ (34) 
dtek-adkQ = ■ ■ ■ (35) 

Boundary conditions for (the incoming modes of) the subsidiary system can be enforced then 
in the standard way. We will consider here for simplicity the 'maximal dissipation' condition, that 
is (we assume that the boundary is on the right): 

G + n^Gk = . (36) 

Now we can use it as a tool for setting up boundary conditions for the energy modes of the main 



evolution system. One can for instance enforce directly (I36p . as in ref. 

We will rather use ()36p as a tool for getting (derivative) boundary conditions for the incoming 
energy mode of the evolution system ([7]- [TT]) . To do this, we can use the evolution equation ^ 



for transforming ()36p into a convenient version of the energy constraint, namely: 

£ = dkV'' -dne + --- (37) 
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We can now use (j37p in order to modify the evolution equation of the incoming energy mode E , 
that is: 



l/a dt E- + dkV'' -dne = a£ + 
The whole process is equivalent to the simple replacement: 

E- ^ E' +a (G^'^'^'') - 9) 
where G*^'^'^'"^ is the solution of the advection equation p6p . 



(38) 



(39) 



The choice a = 1 corresponds to the standard recipe [^-IJ] of ' trading' space normal derivatives 
by time derivatives, in the incoming modes. This implies that the modified mode gets zero propa- 
gation speed along the given direction n. In this case, allowing for (|38p . the time derivative of E~ 
would actually vanish, modulo non-principal terms; this amounts to freezing the incoming modes 
to their initial values (maximal dissipation 'on the right-hand-side'), which is a current practice 
in some Numerical Relativity codes. Note however that constraint preservation requires using the 
right non-principal terms, that can be deduced from the full expression (j39p . 

The choice a = 2 would imply instead that the modified mode gets the same positive speed {+a) 
than the outgoing one E~^ . We show in Appendix B that this choice will lead to a weakly- hyperbolic 
(ill-posed) boundary system. Our results confirm that a = 1 is actually a safe choice Ufl, although 
other values in the interval 1 < a < 2 lead also to a strongly hyperbolic system with non-negative 
speeds for all energy modes (see Appendix B for details). 



Momentum constraint conditions 



The same method can be applied to the momentum constraint modes, although in a less straight- 
forward way. Let us start from the evolution equation (llSp for Zj, and take one extra time deriva- 
tive. We get in this way 

dl Z, + a^dl = ••• (40) 

which, after some cross-derivatives cancellations, leads to the space components of (the principal 
part of) the covariant equation (I3ip . 

A first-order version of (j40p can be obtained again by considering first-order derivatives as 
independent quantities. For the time derivative we will take the obvious choice 

Zi = 1/q dt Zi . (41) 
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The treatment of space derivatives, however, is compHcated by the fact that we are deahng with a 
first-order formulation, so that there are additional ordering constraints to be allowed for. Following 



isl . we will define for further convenience 

Zfc, ^dkZi- ^[fc - ^[fc A] - (1 - C) 7"' Dk]is + (1 + C) 7" A] ks , (42) 



where we have noted for short A = 7'^'' Ars • A closer look to (j42p shows that Z^j is just the space 
derivative of Zj, modulo ordering constraints. In the notation of this paper: 

Zki = -dr ^iki' + % T^fc] + {l + C)[dr f^ikif + ] + • • • (43) 

We can write now (]40p in the first-order form 

dtZi-adkZ'^i = ■■■ (44) 
dt Zki-adk Zi = ■■■ (45) 

which is a symmetric-hyperbolic first-order version of the momentum-constraint evolution system 
(other versions could be obtained by playing with the ordering constraints in a different way). The 
vanishing of the incoming modes of this subsidiary system can be enforced now in the same way 
as for the energy constraint, namely: 

Z, + n'^Zki = . (46) 

This is obviously a 'maximal dissipation' constraint-preserving condition for the subsidiary system, 
which can be used for to get a derivative boundary condition for the main evolution system, as we 
did for the energy modes in the preceding subsection. To be more specific, we can use the evolution 
equation psp for transforming (j46p into a convenient version of the momentum constraint, that is 



M^ = dkU^ + n''Zki + --- (47) 
and use it for modifying the evolution equation of the incoming momentum modes M~ , namely: 

1/q dt Mr + dk X'^m -dnllm = -a Mi + ■ ■ ■ (48) 
which amounts to the following replacement: 

M- ^ M- + a (Zf ^''^ - Zi) , (49) 
where z'f''^^^ is the solution of the advection-like equation (j46p . 
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The choice a = 1 would imply again that the modified modes get zero propagation speeds along 
the normal direction, whereas the choice a = 2 would imply instead that the modified modes get 
the same positive speed (+a) than the outgoing ones . This result requires the extra ordering 
terms in (j42p : this was actually the reason for including them. Note that we can consider different 
values of the coupling parameter for the energy modes (a = a^;), and even for the normal and 
tangent momentum modes {a = , ut , respectively). 

For any value a > 1, the modified modes can be computed consistently from inside. The 
momentum system however is too complicated for a full hyperbolicity analysis, like the one we 
provide in Appendix B for the energy sector. Part of the complication comes from the coupling 
with the non-constraint modes, which require their own boundary conditions. Let us remember 
at this point that the boundary conditions presented in this section are derivative, not algebraic. 
This means that, even in the symmetric hyperbolic cases, proving well-posedness is by no means 
trivial. 

For that reason, we will rather follow the approach of ref. [3]; focusing in the stability of 
small perturbations around smooth solutions, which can be tested numerically. We start in the 
following section, by performing a ' robust stability' test in order to check the numerical stability of 
high-frequency perturbations around the Minkowsky metric. As a full set of boundary conditions is 
required, even in this weak-field test, we supplement our conditions for the constraint-related modes 
with the freezing of the initial values of the incoming non-constraint modes (maximal dissipation 
'on the right-hand-side'). 



IV. NUMERICAL IMPLEMENTATION 



Let us test now the stability and performance of the proposed conditions in the linear weak-field 
regime, by considering a small perturbation of Minkowski space-time, which is generated by taking 
random data both in the extrinsic curvature and in the constraint- violation quantities (G, Zi). 
In this way the initial data violate the energy-momentum constraints, but preserve all ordering 
constraints. The level of the random noise will be of the order 10~^, small enough to make sure that 
we will keep in the linear regime during the whole simulation (Robust Stability test, see ref. 22 1 
for details). 

We will use the standard method of lines [2^ as a finite difference algorithm, so that space and 
time discretization will be treated separately. The time evolution will be dealt with a third-order 
Runge-Kutta algorithm. The time step dt is kept small enough to avoid an excess of numerical 
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dissipation that could distort our results in long runs. 

For space discretization, we will consider a three-dimensional rectangular grid, evenly-spaced 
along every space direction, with a space resolution h = 1/80. We will use there the simplest 
centered, second-order-accurate, discretization scheme. At the points next to the boundary, where 
we can not use the required three-points stencil, we will switch to the standard first-order upwind 
(outgoing) scheme. This combination is the simplest one with the summation-by-parts (SBP) 
property ji^. In this way we expect that the theoretical properties derived from symmetric- 
hyperbolicity will show up in the simulations in a more transparent way. For the same reason, we 
avoid adding extra viscosity terms that could mask the effect of our conditions (at the continuum 
level) with the dissipative effects of the discretization algorithm. Just to make sure, we run also 
with periodic boundary conditions, where the noise level keeps constant: any decrease of the 
constraint-violation level will then be due to the proposed conditions, not to the discretization 
scheme. 

Let us be more specific about the boundary treatment. At boundary points, we use the first- 
order upwind algorithm in order to get a prediction for every dynamical field. Once we have got 
this prediction, we perform the characteristic decomposition along the direction normal to the 
boundary. The predicted values for the outgoing modes, for which the upwind algorithm is known 
to be stable, will be kept (this includes the 'standing' modes, with zero characteristic speed). 
The (unstable) incoming modes will be replaced instead by the values arising from our boundary 
conditions, as described in the preceding section. 

We start with simulations in which the proposed conditions are applied just to the z face, 
whereas we keep periodic boundary conditions along the x and y directions. In this way we can 
detect instabilities which are inherent to the proposed boundary conditions on smooth boundaries 
(no corners), allowing at the same time for some non-trivial dynamics along at least one tangent 
direction. 

We plot in fig. [1] the maximum norm of the energy-constraint-violating quantity G for two 
different choices of the coupling parameter of the energy mode: a^; = 1 , 2 . We can see that, after 
100 crossing times, the case ge = 2 starts showing the effect of the linear modes predicted by 
our hyperbolicity analysis in Appendix B, by departing from the maximal dissipation pattern of 
decay. We plot for comparison the results obtained by applying periodic boundary conditions, so 
we can see how, for the choice a^; = 1, the proposed constraint-preserving conditions are extremely 
effective at 'draining out' energy constraint violations. The rate of decay is actually the same as 
the one obtained by applying maximal dissipation conditions ' on the right-hand-side' also to the 
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FIG. 1: Robust stability test. Time evolution of the maximum norm of O, with just one face open (periodic 
boundaries are implemented along the transverse directions). The fully periodic boundaries result (dashed 
lines) is also included for comparison. We see some growing mode onset in the a^; = 2 case, whereas the 
constraint- preserving — 1 case (continuous line) is very efRcient at reducing the initial noise level. 

energy modes, as expected from the analysis given in the previous section. In what follows, we will 
fix a£; = 1 for this coupling parameter. 

We plot in fig. [2] both the maximum norm of the longitudinal (left panel) and transverse com- 
ponents (right panel) of the momentum-constraint-violating vector Zj for the choice a^r = ay = 1 





0.1 1 10 100 0.1 1 10 100 



Crossing times Crossing times 

FIG. 2: Robust stability test. Time evolution of the maximum norm of the longitudinal and transverse Zi 
components (left and right panels, respectively). In both cases, the periodic boundaries results (dashed lines) 
are included for comparison. The initial noise in the momentum constraint gets reduced very efficiently in 
both the C = — 1 and the C — cases, although there is a slight difference, more visible in the longitudinal 
case (left panel). 
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of the coupling parameter of the momentum modes. We include again for comparison the results 
obtained by applying periodic boundary conditions, so we can see how the proposed constraint- 
preserving conditions are very effective at 'draining out' energy constraint violations. The Zy plots 
are slightly, sensitive to the ordering parameter = 0, — 1. In the ^ = — 1 case, the rate of decay 
is actually the same as the one obtained by applying instead maximal dissipation conditions ' on 
the right-hand. -side' for the momentum modes. The results are qualitatively the same for other 
components of Zj and for other parameter combinations otv, = 1, 2 . 
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FIG. 3: Robust stability test. Time evolution of the maximum norm of the constraint-violating quantities 
(left panel), and Zy (right panel). The proposed boundary conditions are applied here to all faces, including 
corners and vertices. Some amount of numerical dissipation has been added, so that the periodic boundaries 
plots (dashed lines) get a visible negative slope. The choice = 1 for the energy modes is still clearly 
stable (left panel). The choice C = ~1 for the momentum modes (right panel) shows a growing mode onset. 
For comparison, a plot with the maximal dissipation results is also included in the right panel (bottom line). 



In order to perform a full test for cartesian-like grids, including corner and vertex points, we 
will repeat the same simulations, but this time with the proposed boundary conditions applied to 
all faces, not iust to the z ones. A standard treatment of corner points 'a la Olsson' ITil, like the 
one presented in previous works [Slilla]) results into numerical instability issues. A simple cure is to 
add some extra dissipation at the interior points, at the price of masking the theoretical results, et 
the continuum level, with the numerical viscosity effects, as shown in fig. [31 We can see there that 
opening all faces makes the effects to appear much faster. The expected instability of the a^; = 2 
choice of the energy coupling parameter, which was just an onset in fig. [H shows up manifestly 
here (left panel). Also, a growing mode onset is clearly visible for the choice C, = —\ oi the 
momentum-constraint coupling parameter (right panel). The case C = looks stable, although no 
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strong conclusion can be drawn because of the added numerical dissipation. Maximal dissipation 

results are also shown for comparison in the right panel (bottom line). 

• , ^ 



^ . 1 t 

1 1 1 1 1 ^ 

I 1 1 1 1 i 

! L [ 1 J i 

FIG. 4: Stencil for first-order prediction at boundary points (black dots). Interior points belong to the 
shaded zone. The stencil for two boundary points is represented by thick lines (one space direction has been 
suppressed for clarity). Note that no corner points are required, as the tangent derivatives are not computed 
at the boundary, but at the neighbor layer. 

We will present here an alternative numerical treatment. At boundary points, tangent deriva- 
tives are computed at the next-to-last layer. The corresponding stencil is shown in fig. |4l In this 
way the corner points are not required. This avoids the reported code stability issues, even with- 
out adding extra numerical dissipation terms. Note that transverse derivatives are still computed 
using the standard three-point SBP algorithm, like in the smooth boundaries case. As every space 
derivative can be considered separately (we are dealing with a first-order system) the SBP property 
should still follow for our modified scheme. The price for the shift of the transverse derivatives 
to the next-to-last layer is getting just first-order accuracy at the boundary, but the longitudinal 
derivatives there were yet only first-order accurate anyway. 

This discretization variant allows getting stable results, at least for the value C = of the 
ordering parameter. We plot in fig. [5] the maximum norm of the constraint-violation quantities 
(0, Zi). We can see there that removing the extra numerical dissipation makes the features more 
transparent. The instability of the a^; = 2 choice of the energy coupling parameter, appears now 
instantly. The downfall rate in the stable case a^; = 1, increased as the constraint violations are 
drained out in all three directions now, can be seen in a more unambiguous way. Concerning the 
momentum constraint (right panel), the standard ^ = — 1 ordering shows now clearly its unstable 
behavior, which was masked by the added dissipation in the standard treatment (see fig. [3|). The 
centered ordering choice C = recovers instead the manifest stable behavior shown in single-face 
simulations (see fig. [2]), close to the maximal dissipation case (left panel, bottom line). 
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FIG. 5: Robust stability test. Time evolution of the maximum norm of the constraint-violating quantities O 
(left panel), and Zy (right panel). The proposed boundary conditions are applied to all faces. Corner points 
are avoided in the way shown in fig. |4l No extra numerical dissipation has been added, so that the periodic 
boundaries plots (dashed lines) keep flat. The absence of extra dissipation clarifies the features shown in 
the previous figure. 

Our results show the numerical stability of the proposed boundary conditions in the linear regime 
for suitable combinations of the coupling and/or ordering parameters. The proposed boundary 
conditions produce instead a very effective decreasing of (the cumulated effect of) energy and 
momentum constraint violations which compares with the one obtained by applying maximal 
dissipation boundary conditions to (the right-hand-side of) the constraint related modes. 

V. GOWDY WAVES AS A STRONG FIELD TEST 



Although the results of the preceding are encouraging, let us remark that we were just testing 
the linear regime around Minkowsky spacetime. This is not enough, as high-frequency instabilities 



can appear in generic, strong field, situations (see for instance ref. 



). In order to test the strong 



20|, which describes a space-time containing 



field regime, we will consider the Gowdy solution 
plane polarized gravitational waves. This is one of the test cases that is used in numerical code 
cross-comparison with periodic boundary conditions {22 , 



23( . One of the advantages is that it allows 



for periodic and/or reflecting boundary conditions, which can be applied to the modes which are 
not in the energy- momentum constraint sector. A first proposal for this selective testing of the 
constraint-related modes has been presented recently 211]. 
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The Gowdy line element can be written as 

ds' = e«/2 i-dt' + dz^) + t (e^ dx^ + e-^ dy^) (50) 

n 

where the quantities Q and P are functions of t and z only and periodic in z, that is ^22i] . 

P = Jo{2TTt) cos(27rz) (51) 
Q = 7rJo(27r)Ji(27r) - 27rtJo(27rt)Ji(27rt)cos^(27rz) 

+ 27r2t2[Jo2(27rt) + Ji2(27rt)- Jo2(27r)- Ji2(27r)] (52) 

so that the lapse function a = t~^/^ e'^/^ is constant in space at any time to at which Jo(27rio) 
vanishes. 

Let us now perform the following time coordinate transformation 

t = toe-^'^\ (53) 

so that the expanding line element (|50p is seen in the new time coordinate r as collapsing towards 
the t = singularity, which is approached only in the limit r — > oo. This 'singularity avoidance' 
re r coordinate follows from the fact that the resulting slicing by r = const surfaces 



property of t 

is harmonic 25l]. We will launch our simulations in normal coordinates, starting with a constant 
lapse aQ = 1 Sit T = Q {t = to). 

The discretization is performed like in the preceding section, but with a space resolution h = 
1/100. Allowing for the fact that the only non-trivial space dependence in the metric is through 
cos(27rz), the numerical grid is fitted to the range < z < 1. In this way, the exact solution admits 
either periodic or reflection boundary conditions. We can use these exact boundary conditions as 
a comparison with the constraint-preserving ones that we are going to test. As the Gowdy metric 
components depend on just one coordinate, we will apply here the proposed constraint-preserving 
conditions only to the z faces, keeping periodic boundary conditions along the transverse directions. 
Also, like in the preceding section, we show the results for the aE = = clt = ^ coupling 
parameters combination, although other choices of ajv, ax = 1; 2 lead to similar results. 

We start with a simple convergence test. As we know the exact solution (jSOh . we can directly 
compute the relative error of every simulation. Then, only two different resolutions are required for 
checking convergence. We will take h = 1/50 and h = 1/100 for our Gowdy wave simulations. We 
plot in fig. [6] the energy-constraint- violation quantity at some early time t = 10 (left panel). We 
see the expected second-order accuracy at the interior points which are yet causally disconnected 
with the boundaries. We see just first-order accuracy at the boundary, plus a smooth transition 
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FIG. 6: Convergence test. The constraint-violating quantity 9 is plotted at i = 10 for two different 
resolutions (left panel) . We see second-order convergence in the interior region, decreasing up to first-order 
at points causally connected to the boundary, as expected from our SBP algorithm. We plot in the right 
panel the relative error of the gzz metric component, evolved up to t = 250. The boundary-induced accuracy 
reduction is not even visible yet. 



zone. This accuracy reduction at boundaries is inherent to simple SBP algorithms, which require a 
lower-order discretization at boundary points [17]. One could keep instead the accuracy level of the 
interior points by using more accurate predictions for boundary values, but at the price of loosing 
the SBP property. In our test case, however, this issue is not affecting the metric components, 
even at a much later time t = 250, as we can see in the right panel of fig. [H 
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FIG. 7: Gowdy waves test. The 9 and profiles are plotted as indicators of the accumulated error 
due to energy-momentum constraint violations. Reflection boundaries results are also plot for comparison 
(continuous lines). Dotted lines correspond to the proposed boundary conditions, whereas dashed lines 
correspond to the same conditions with the extra damping terms discussed in this section, with rj = 0.1. 
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We show in fig. [7] the results of the h = 1/100 resolution simulation for the boundary profiles 
of Q and Zi, indicating the accumulated amount of energy and momentum constraint violations 
(left and right panels, respectively). We apply in both cases the proposed boundary conditions at 
the z faces to the constraint-related modes, while keeping exact (reflection) boundary conditions 
for the other modes. We can see that the constraint-preserving conditions result, in this strong- 
field test, into an accumulated amount of constraint violations (dotted lines) that is similar or 
even slightly better than the one produced by the interior points treatment, which can be seen 
in simulations with (exact) reflection boundaries for all faces (continuous lines). Note that the 
reflection conditions anchor to zero at the boundary points, which is always more accurate in 
this test, although not very useful in more realistic cases. 

These results confirm that the proposed boundary conditions are indeed constraint-preserving, 
in the sense that their contribution to energy and momentum constraint violations keeps within the 
limits of the truncation error of the discretization algorithms, even in this strong field scenario. This 
good behavior can be further improved by introducing constraint-damping terms in the evolution 
equations for the boundary quantities (f36l H6]l that is 

e + n^Gfc = -r/ e (54) 
Zi + n^^Zki = -vZi. (55) 

The resulting values can then be used in the replacements (139p and (149p . respectively. We have in- 
cluded the corresponding results in fig. [7] (dashed lines). The amount of both energy and constraint 
violations becomes even lower than the one for the (exact) reflection boundaries simulations even 
for a small value rj = 0.1 of the damping parameter. The effect is specially visible in the energy 
constraint case (left panel). 



VI. CONCLUSIONS AND OUTLOOK 

The work presented in this paper revises and improves the previous results for the Z4 case [sl-fl^ 
in many different ways. 

On the theoretical side, we have proposed a new symmetrizer, which extends the parametric 
domain for symmetric hyperbolicity from the single value C = —1 to the interval — 1 < C ^ 0- 
We have identified in the process a new basis for the dynamical field space (l2][5]) which allows 
a clear-cut separation between the constraint-related modes and the remaining ones. Regarding 
the boundary treatment, we have also generalized the way in which boundary conditions can by 
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used for modifying the incoming modes, by introducing a new parameter a which, at least for 
the momentum constraint modes, can depart from the standard value a = 1 without affecting the 
stability of the results. 

On the numerical side, the use of the new basis definitely improves the stability of the previous 
Z4 results. In the single face case, where we use periodic boundary conditions along transverse 
directions, we see that the linear modes previously reported in the robust stability test P, [3] for 
the symmetric ordering case ((" = 0) no longer show up. Moreover, we have devised a simple finite- 
differences stencil for the prediction step at the boundaries which avoids the corner and vertex 
points even in cartesian-like grids, providing an interesting alternative to the standard (Olsson) 
corners treatment. 

The proposed boundary conditions have been also tested in a strong field scenario, the Gowdy 
waves metric so that the effect of non-trivial metric coefficients can be seen in the simulation 
results. The convergence test in this non-linear regime provides strong evidence of numerical stabil- 
ity for some suitable parameter combinations. Our simulations actually confirm that the proposed 
boundary conditions are constraint-preserving: the accumulated amount of energy-momentum con- 
straint violations is similar or even better than the one generated by either periodic or reflection 
conditions, which are exact in the Gowdy waves case. 

Now it remains the question of how these interesting results can be extended to other 3-1-1 evo- 
lution formalisms and/or gauge conditions. Let us remember that all our symmetric hyperbolicity 
results apply as usual just to the harmonic slicing, not to the 'l-|-log' class of slicing conditions 
which are customary in BSSN black-hole simulations. There is no problem, however, in extending 
the proposed boundary conditions to this case: in our new basis the gauge sector is clearly sepa- 
rated from the constraint-related one, so that one can keep using the replacements l\39\ H9|) even in 
this non-harmonic case. The shift, however, introduces new couplings and would require a detailed 
case-by-case investigation: even the strong hyperbolicity of the system can depend on the specific 
choice of shift condition. 

Concerning the extension from the Z4 to the BSSN formalism, the momentum constraint treat- 
ment can be derived from the simple condition [l^ 

fi = -7ifc7'^j + 2Z* (56) 

which relates the additional BSSN quantity Fj with the space vector Zj. The replacement ()49p 
can then be used for getting a suitable boundary condition in this context. The case of the energy 
constraint is more challenging, as the BSSN formalism does not contain any supplementary quantity 
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analogous to Q. One could follow, however, the line recently proposed in [l^: a slight modification 
of the original BSSN equations allows to include the new quantity 0, so that the correspondence 
with the Z4 formalism is complete. The replacement (I39p can then be used directly in such context. 

A major challenge is posed by the fact that most BSSN implementations are of second order in 
space. This has some advantages in this context, as the ordering constraints do not show up and 
this removes the main source of ambiguities in the constraint-violations evolution system. As a 
result, the boundary conditions ([36l H6]) become simply advection equations so that we can expect 
a more effective constraint-violation 'draining' rate. The problem, however, is that second-order 
implementations do not have the algebraic characteristic decomposition which is crucial in the first- 
order ones. The boundaries treatment takes quite different approaches in second-order formalisms, 
although the evolution equations for the constraint-related quantities Q , Zi are still of first order 
in the Z4 case, even at the continuum level, and this suggests that the results presented here can 
be still helpful in this case. We are currently working in this direction. 
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Appendix A: Positivity of the energy estimate 

We have derived in section II an generalized ' energy estimate' for the Z4 system, namely: 

s = e^ + VkV'^ + u'mij + ~^^'^jikij + (i + c){z^Zk - jJ^'^hjk) + 2 c z^w^ , (a.i) 

where we noted 

l^kij = fJ-kij - Wk Jij ■ (A. 2) 

In order to check the positivity of (jA.ip . let us consider the decomposition of the three-index 
tensor /i^jj into its symmetric and antisymmetric parts, that is 

f^kij = fi-ikij) + f^tij ■ (^-3) 
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Allowing for the identities, 

,-.(") _n ,-.(") _ _i r.W (A A) 

the rank-three terms contribution to S can the be written as 

s=-c A^'^^') + ^^ /i^") + • • • (A.5) 

(the dots stand for lower-rank components). It follows that a necessary condition for positivity is 
> C > -3. 

We can now rewrite (|A.ip as 

5 = 02 + y^yfc + n^i n,, - c /ifc., A''^' + ^ (1 + /ig /i^"^ + (1 + C) ^'^fc + 2 C z,,w\ (A.6) 

Allowing for ()A.2p . which implies in turn 

/i\, = -Z,-Ty,, (A.7) 
we see that we can rewrite again ()A.6P as 

5 = + + tfm,, - c Xmj 'x"'' + ^ (1 + C) /ig /^^"^ + (i + C) ^'^^ , (a.s) 

where 

Xkij = Xkij lc=-i = fj-kij + lk{iWj) ■ (A-9) 

It follows from the final expression (jA.Sp that the energy estimate is positive definite in the whole 
interval 

> C > - 1- (A.IO) 
Note that for = —1, that is A = A, we recover the estimate given in ref. [l^. 

Appendix B: Strong hyperbolicity of the energy modes 

We can analyze the hyperbolicity of the boundary evolution system, by considering the char- 
acteristic matrix along a generic oblique direction r, which is related to the normal direction n 
by 

r = n coscp + s sinip , (B.l) 
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where we have taken 

= = 1 n • s = . (B.2) 

The strong hyperbohcity requirement amounts to demand that the characteristic matrix is fully 
diagonalizable and has real eigenvalues (propagation speeds) for any value of the angle ^p. 

In order to compute the characteristic matrix, we will consider the standard form of (the prin- 
cipal part of) the evolution system as follows 

dt u + a5rF''(u) = ••• , (B.3) 

where u stands for the array of dynamical fields and is the array of fluxes along the direction 
r. We will restrict ourselves here to the Energy-modes subsystem, which consists in the fields 

u = (£+ , i^- , , T/p ) (B.4) 

the index p meaning here a projection along the direction orthogonal both to n and s. 

It is clear that the two components Vp are eigenvectors of the characteristic matrix with zero 
propagation speed. The non-trivial fiuxes are then: 

F''(S+) = cosLp + sin^ Vs (B.5) 
F''{E-) = {a - I) cosip E- + {I - a) simp Vs (B.6) 
F'{Vs) = \ sin^ {E+ + E~), (B.7) 

where we have allowed for the modified evolution equation (|38p . We can see that the one of the 
characteristic speeds is zero and the other two are be given by the solutions of 

{v — a cosip) (f — (o — 1) a cosip) = (1 — -) sin^^p , (B.8) 

which has real distinct solutions for a < 2 . The degenerate case a = 2 is not diagonalizable. It 
follows that the boundary evolution subsystem given by the above fluxes is strongly hyperbolic for 
a < 2 and weakly hyperbolic for a = 2 . 
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